package fem2.element.crack_sda2;

import inf.jlinalg.Array1DMatrix;
import inf.jlinalg.Array2DMatrix;
import inf.jlinalg.ArrayVector;
import inf.jlinalg.BLAM;
import inf.jlinalg.IMatrix;
import inf.jlinalg.IMatrixReadOnly;
import inf.jlinalg.IVector;
import inf.jlinalg.IVectorReadOnly;
import inf.math.LinAlg;

import java.util.ArrayList;

import math2.LinAlgUtilities;
import fem2.Debugger;
import fem2.IntegrationManager;
import fem2.MaterialModel;
import fem2.MeshPart;
import fem2.element.StructuralElement;
import fem2.enu.EchoLevelType;
import fem2.enu.MeshPartType;
import fem2.material.cohesive_law2.CohesiveLaw;

/**
 * Strong discontinuity approach element
 * 
 * @author hbui
 * 
 */
public abstract class SDAElement extends StructuralElement {

	protected ArrayList<double[]> crackBoundary;

	protected double[] uHatLocal; // buffer uHatLocal

	private boolean isCracked;

	private int p_crack;

	private MeshPartType crackMeshPartType = MeshPartType.LINE2;

	private CohesiveLaw cl;

	protected IVector alpha; // enhanced strain mode

	public SDAElement(MeshPart mp, MaterialModel mm, CohesiveLaw cl) {
		super(mp, mm);
		setCohesiveLaw(cl);
	}

	/**
	 * set the crack mesh part type. Crack mesh part type can be 2-node line or
	 * face
	 * 
	 * @param type
	 */
	public void setCrackMeshPartType(MeshPartType type) {
		this.crackMeshPartType = type;
	}

	/**
	 * set the cohesive law (i.e crack material model)
	 * 
	 * @param cl
	 */
	public void setCohesiveLaw(CohesiveLaw cl) {
		this.cl = cl;
		if (this.cl != null) {
			this.cl.registerElement(this);
		}
	}

	/**
	 * get the cohesive law associated with the element
	 * 
	 * @return
	 */
	public CohesiveLaw getCohesiveLaw() {
		return this.cl;
	}

	/**
	 * get the mesh part type of the crack
	 * 
	 * @return
	 */
	public MeshPartType getCrackMeshPartType() {
		return crackMeshPartType;
	}

	/**
	 * 
	 * @return true if element is cracked; false otherwise
	 */
	public boolean isCracked() {
		return this.isCracked;
	}

	/**
	 * get the crack opening parameters associated with the element
	 * 
	 * @return
	 */
	public IVector getCrackOpeningParameter() {
		return alpha;
	}

	/**
	 * get lists of jump boundary points
	 * 
	 * @return
	 */
	public ArrayList<double[]> getX() {
		return crackBoundary;
	}

	/**
	 * 
	 * @return the number of crack opening parameters associated with the
	 *         element
	 */
	public int getNumberOfCrackOpeningParameters() {
		return getCrackOpeningParameter().getSize();
	}

	/**
	 * set integration order on crack surface
	 * 
	 * @param p
	 */
	public void setIntegrationOrderOnCrack(int p) {
		this.p_crack = p;
	}

	/**
	 * 
	 * @return number of Gauss point on the crack
	 */
	public int getNumGaussPointOnCrack() {
		return IntegrationManager.getInstance().getNumIntegrationPoints(getCrackMeshPartType(),
				p_crack);
	}

	/**
	 * 
	 * @param k
	 * @return Gauss point on crack
	 */
	public double[] getGaussPointOnCrack(int k) {
		return IntegrationManager.getInstance().getIntegrationPoint(getCrackMeshPartType(),
				p_crack, k);
	}

	/**
	 * 
	 * @param k
	 * @return Gauss weight on crack
	 */
	public double getGaussWeightOnCrack(int k) {
		return IntegrationManager.getInstance().getIntegrationWeight(getCrackMeshPartType(),
				p_crack, k);
	}

	/**
	 * check crack onset criteria
	 * 
	 * @return true if crack is initiated; otherwise is false
	 */
	protected abstract boolean checkOnsetCriteria();

	/**
	 * compute determinant of Jacobian on crack surface
	 * 
	 * @param s
	 * @return
	 */
	protected abstract double computeDeterminantOfCrackJacobian(double... s);

	/**
	 * Compute enhanced strain operator using in compute sigma. sigma = Bu +
	 * G(xi) * alpha
	 * 
	 * @param mp
	 * @param g
	 * @param xi
	 * @return the enhanced strain operator
	 */
	protected abstract void computeEnhancedStrainOperator(MeshPart mp, IMatrix g, double... xi);

	/**
	 * Compute the projection operator for sigma in computing the local
	 * equilibrium. rlocal = int(G(x)^T * sigma, \Omega) + int(t, \Gamma)
	 * 
	 * @param mp
	 * @param g
	 * @param xi
	 * @return the projection operator
	 */
	protected abstract void computeProjectionOperator(MeshPart mp, IMatrix g, double... xi);

	/**
	 * compute crack traction on crack surface
	 * 
	 * @param k
	 *            index of intergration point on crack
	 * @return
	 */
	protected abstract IVector computeCrackTraction(int k);

	/**
	 * compute crack stiffness
	 * 
	 * @param k
	 *            index of intergration point on crack
	 * @return
	 */
	protected abstract IMatrix computeCrackStiffness(int k);

	@Override
	public IVector computeStrain(double... xi) {
		if (!isCracked()) {
			return super.computeRegularStrain(xi);
		} else {
			MeshPart mp = getMeshPart();
			int dim = mp.getDimension();
			int ntrs = dim * (dim + 1) / 2;
			int na = getNumberOfCrackOpeningParameters();
			IMatrix G_enhanced = Array1DMatrix.createColumnMajor(ntrs, na);
			computeEnhancedStrainOperator(mp, G_enhanced, xi);
			IVector alpha = getCrackOpeningParameter();
			IVector e = super.computeRegularStrain(xi);
			BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, G_enhanced, alpha, 1.0, e);
			return e;
		}
	}

	@Override
	public IMatrixReadOnly computeKte(IMatrix ke) {
		if (!isCracked()) {
			ke = (IMatrix) super.computeKte(ke);
			if (getEchoLevel().contain(EchoLevelType.AUX2)) {
				Debugger.watch("Kte = ", ke);
			}
			return ke;
		} else {
			MeshPart mp = getMeshPart();
			int n = mp.getBasis().getDimension();
			int ng = getNumIntegrationPoints();
			int dim = mp.getDimension();
			int ntrs = dim * (dim + 1) / 2;
			int na = getNumberOfCrackOpeningParameters();
			double th = getMaterialModel().getSectionProperty();
			Array1DMatrix B = Array1DMatrix.createColumnMajor(ntrs, dim * n);
			IMatrix C = null;
			IMatrix G_enhanced = Array1DMatrix.createColumnMajor(ntrs, na);
			IMatrix G_project = Array1DMatrix.createColumnMajor(ntrs, na);

			IMatrix Kbb = ke;
			IMatrix Kbg = new Array2DMatrix(dim * n, na);
			Array1DMatrix Kgg = Array1DMatrix.createColumnMajor(na, na);
			IMatrix Kgb = new Array2DMatrix(na, dim * n);

			for (int k = 0; k < ng; k++) {
				double w = getIntegrationWeight(k);
				double[] xi = getIntegrationPoint(k);
				double[][] jac = mp.jacXAt(xi);
				double detJ = LinAlgUtilities.determinantNxN(jac);
				double coeff = w * detJ * th;

				computeB(mp, B, xi);
				C = getMaterialModel().getConstitutiveMatrix(this, k);
				computeEnhancedStrainOperator(mp, G_enhanced, xi);
				computeProjectionOperator(mp, G_project, xi);

				LinAlgUtilities.btCB(B, C, coeff, Kbb);

				LinAlgUtilities.atCB(B, G_enhanced, C, coeff, Kbg);

				LinAlgUtilities.atCB(G_project, B, C, coeff, Kgb);

				LinAlgUtilities.atCB(G_project, G_enhanced, C, coeff, Kgg);
			}

			if (getEchoLevel().contain(EchoLevelType.AUX3)) {
				Debugger.watch("Kgg_1 = ", Kgg);
			}

			int ngc = getNumGaussPointOnCrack();
			for (int k = 0; k < ngc; k++) {
				double w = getGaussWeightOnCrack(k);
				double[] s = getGaussPointOnCrack(k);
				double detJ = computeDeterminantOfCrackJacobian(s);
				double coeff = w * detJ * th;
				IMatrix Cgamma = computeCrackStiffness(k);
				BLAM.add(coeff, Cgamma, Kgg);

				if (getEchoLevel().contain(EchoLevelType.AUX3)) {
					Debugger.watch("Cgamma = ", Cgamma);
				}

			}

			if (getEchoLevel().contain(EchoLevelType.AUX3)) {
				Debugger.watch("last C = ", C);
				Debugger.watch("last B = ", B);
				Debugger.watch("Kbg = ", Kbg);
				Debugger.watch("Kgb = ", Kgb);
				Debugger.watch("Kgg = ", Kgg);
			}

			LinAlgUtilities.backSolve(Kgg, Kgb);
			BLAM.multiply(-1.0, BLAM.NO_TRANSPOSE, Kbg, BLAM.NO_TRANSPOSE, Kgb, 1.0, Kbb);

			if (getEchoLevel().contain(EchoLevelType.AUX2)) {
				Debugger.watch("Kte = ", ke);
			}

			return Kbb;
		}
	}

	@Override
	public IVectorReadOnly computeRie(IVector rie) {
		if (!isCracked()) {
			rie = (IVector) super.computeRie(rie);
			if (getEchoLevel().contain(EchoLevelType.AUX2)) {
				Debugger.watch("rie = ", rie);
			}
			return rie;
		} else {
			MeshPart mp = getMeshPart();
			int n = mp.getBasis().getDimension();
			int ng = getNumIntegrationPoints();
			int dim = mp.getDimension();
			int ntrs = dim * (dim + 1) / 2;
			int na = getNumberOfCrackOpeningParameters();
			double th = getMaterialModel().getSectionProperty();
			Array1DMatrix B = Array1DMatrix.createColumnMajor(ntrs, dim * n);
			IMatrix G_enhanced = Array1DMatrix.createColumnMajor(ntrs, na);
			IMatrix G_project = Array1DMatrix.createColumnMajor(ntrs, na);

			IMatrix Kbg = new Array2DMatrix(dim * n, na);
			Array1DMatrix Kgg = Array1DMatrix.createColumnMajor(na, na);
			ArrayVector riLocal = new ArrayVector(na);

			for (int k = 0; k < ng; k++) {
				double w = getIntegrationWeight(k);
				double[] xi = getIntegrationPoint(k);
				double[][] jac = mp.jacXAt(xi);
				double detJ = LinAlgUtilities.determinantNxN(jac);
				double coeff = w * detJ * th;

				computeB(mp, B, xi);
				IVector s = getStress(k);
				IMatrix C = getMaterialModel().getConstitutiveMatrix(this, k);
				computeEnhancedStrainOperator(mp, G_enhanced, xi);
				computeProjectionOperator(mp, G_project, xi);

				BLAM.multiply(coeff, BLAM.TRANSPOSE, B, s, 1.0, rie);

				LinAlgUtilities.atCB(B, G_enhanced, C, coeff, Kbg);

				LinAlgUtilities.atCB(G_project, G_enhanced, C, coeff, Kgg);

				BLAM.multiply(coeff, BLAM.TRANSPOSE, G_project, s, 1.0, riLocal);
			}

			int ngc = getNumGaussPointOnCrack();
			for (int k = 0; k < ngc; k++) {
				double w = getGaussWeightOnCrack(k);
				double[] s = getGaussPointOnCrack(k);
				double detJ = computeDeterminantOfCrackJacobian(s);
				double coeff = w * detJ * th;
				IMatrix Cgamma = computeCrackStiffness(k);
				BLAM.add(coeff, Cgamma, Kgg);
				IVector t = computeCrackTraction(k);
				BLAM.add(coeff, t, riLocal);
			}

			LinAlgUtilities.backSolve(Kgg, riLocal.getValues());
			BLAM.multiply(-1.0, BLAM.NO_TRANSPOSE, Kbg, riLocal, 1.0, rie);

			if (getEchoLevel().contain(EchoLevelType.AUX2)) {
				Debugger.watch("rie = ", rie);
			}
			return rie;
		}
	}

	@Override
	public void Initialize() {
		super.Initialize();
		cl.Initialize(this);
		uHatLocal = new double[getMeshPart().countDofs()];
	}

	@Override
	public void InitializeSolutionStep() {
		super.InitializeSolutionStep();
		cl.InitializeSolutionStep(this);
		/*
		 * save ue
		 */
		LinAlg.copy(uHatLocal.length, getMeshPart().getUHatLocal(), uHatLocal);
	}

	@Override
	public void InitalizeNonlinearIteration() {
		super.InitalizeNonlinearIteration();
		cl.InitalizeNonlinearIteration(this);
	}

	@Override
	public void FinalizeNonlinearIteration() {
		if (isCracked()) {
			/*
			 * Note: this type of element use static consensation to condense
			 * the jump enhanced mode. Therefore, we must make sure the riLocal
			 * and other condensed matrices are the same in computeRi &
			 * computeKte. Storing matrices/vectors is a more credible solution
			 * but it's not memory-efficient for big model. Therefore
			 * recomputation is needed but the consistency must be taken care.
			 */
			MeshPart mp = getMeshPart();
			int n = mp.getBasis().getDimension();
			int ng = getNumIntegrationPoints();
			int dim = mp.getDimension();
			int ntrs = dim * (dim + 1) / 2;
			int na = getNumberOfCrackOpeningParameters();
			double th = getMaterialModel().getSectionProperty();
			Array1DMatrix B = Array1DMatrix.createColumnMajor(ntrs, dim * n);
			IMatrix G_enhanced = Array1DMatrix.createColumnMajor(ntrs, na);
			IMatrix G_project = Array1DMatrix.createColumnMajor(ntrs, na);

			ArrayVector riLocal = new ArrayVector(na);
			Array1DMatrix Kgg = Array1DMatrix.createColumnMajor(na, na);
			IMatrix Kgb = new Array2DMatrix(na, dim * n);

			for (int k = 0; k < ng; k++) {
				double w = getIntegrationWeight(k);
				double[] xi = getIntegrationPoint(k);
				double[][] jac = mp.jacXAt(xi);
				double detJ = LinAlgUtilities.determinantNxN(jac);
				double coeff = w * detJ * th;

				computeB(mp, B, xi);
				IVector sigma = getStress(k);
				IMatrix C = getMaterialModel().getConstitutiveMatrix(this, k);
				computeEnhancedStrainOperator(mp, G_enhanced, xi);
				computeProjectionOperator(mp, G_project, xi);

				LinAlgUtilities.atCB(G_project, G_enhanced, C, coeff, Kgg);

				LinAlgUtilities.atCB(G_project, B, C, coeff, Kgb);

				BLAM.multiply(coeff, BLAM.TRANSPOSE, G_project, sigma, 1.0, riLocal);

				if (getEchoLevel().contain(EchoLevelType.AUX3)) {
					Debugger.watch("G_project = ", G_project);
					Debugger.watch("G_enhanced = ", G_enhanced);
				}
				if (getEchoLevel().contain(EchoLevelType.AUX4)) {
					Debugger.watch("sigma = ", sigma);
				}
			}

			int ngc = getNumGaussPointOnCrack();
			for (int k = 0; k < ngc; k++) {
				double w = getGaussWeightOnCrack(k);
				double[] s = getGaussPointOnCrack(k);
				double detJ = computeDeterminantOfCrackJacobian(s);
				double coeff = w * detJ * th;
				IMatrix Cgamma = computeCrackStiffness(k);
				BLAM.add(coeff, Cgamma, Kgg);
				IVector t = computeCrackTraction(k);

				if (getEchoLevel().contain(EchoLevelType.AUX4)) {
					Debugger.watch("traction = ", t);
				}

				BLAM.add(coeff, t, riLocal);
			}

			ArrayVector due = new ArrayVector(mp.getUHatLocal().clone());
			LinAlg.add(uHatLocal.length, -1.0, uHatLocal, 1.0, due.getValues());

			// if (getEchoLevel().contain(EchoLevelType.AUX3)) {
			// Debugger.watch("Kgg = ", Kgg);
			// Debugger.watch("Kgb = ", Kgb);
			// }
			if (getEchoLevel().contain(EchoLevelType.AUX4)) {
				Debugger.watch("riLocal = ", riLocal);
				Debugger.watch("Xi = ", getCrackOpeningParameter());
				Debugger.watch("ue = ", mp.getUHatLocal());
				Debugger.watch("ue_old = ", uHatLocal);
				Debugger.watch("due = ", due);
			}

			BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, Kgb, due, 1.0, riLocal);
			LinAlgUtilities.backSolve(Kgg, riLocal.getValues());

			/*
			 * update crack opening parameter
			 */
			IVector alpha = getCrackOpeningParameter();
			BLAM.add(-1.0, riLocal, alpha);

			if (getEchoLevel().contain(EchoLevelType.AUX4)) {
				Debugger.watch("dalpha = -", riLocal);
			}

			/*
			 * update cohesive law
			 */
			cl.FinalizeNonlinearIteration(this);
		}

		/*
		 * update material
		 */
		super.FinalizeNonlinearIteration();

		/*
		 * save ue
		 */
		LinAlg.copy(uHatLocal.length, getMeshPart().getUHatLocal(), uHatLocal);
	}

	@Override
	public void FinalizeSolutionStep() {
		if (!isCracked) {
			isCracked = checkOnsetCriteria();
			if (isCracked) {
				Debugger.warn("Crack initiated at element " + getId());
			}
		}
		/*
		 * output & debugging
		 */
		if (isCracked()) {
			// /*
			// * check equilibrium condition
			// */
			// double norm_rilocal = BLAM.norm(BLAM.NORM_TWO, rilocal);
			// if (norm_rilocal > MyLinAlg.tol6) {
			// Debugger.error("crack equilibrium condition is not met at " +
			// getElemType() + " "
			// + getMeshPart().getId() + ", norm(renh) = " + norm_rilocal);
			// }
			/*
			 * print out results after step complete
			 */
			if (getEchoLevel().contain(EchoLevelType.OUTPUT)) {
				Debugger.watch("crack opening parameters: ", getCrackOpeningParameter());
			}
			if (getEchoLevel().contain(EchoLevelType.AUX3)) {
				for (int k = 0; k < getNumIntegrationPoints(); k++) {
					Debugger.watch("sigma " + k + " =", getStress(k));
				}
			}

			/*
			 * only update when element is cracked
			 */
			cl.FinalizeSolutionStep(this);
		}

		super.FinalizeSolutionStep();
	}

	@Override
	public void CalculateOutputData() {
		/*
		 * TODO
		 */
	}
}
